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A  review  is  presented  for  two-phase  modeling  approaches  to  study  various  transport  processes  and 
reactions  in  polymer  electrolyte  membrane  (PEM)  fuel  cells  along  with  some  experimental  work.  It  has 
been  noted  that  water  management  is  still  one  of  the  least  accurate  modeled  phenomena.  The  lackness 
in  complete  descriptive  models  for  water  management  inside  PEM  fuel  cells  can  be  attributed  to  the 
complexity  of  the  phenomena,  lack  of  empirical  or  measured  data  and  non-availability  of  apt  governing 
equations. 

Another  discrepancy  found  in  present  models  is  the  proper  validation  of  the  numerical  work  as  it 
has  been  observed  that  mere  comparison  with  V-I  curve  can  sometimes  lead  to  misguided  conclusions. 
Additionally,  keeping  in  mind  the  multi-scale  nature  of  a  PEM  fuel  cell,  application  of  the  Lattice  Boltz¬ 
mann  (LB)  method  has  also  been  reviewed  in  this  work  and  it  was  noticed  that  LB  methods  offer  bright 
perspective  at  meso-scale  by  incorporating  details  of  local  structure.  Furthermore,  a  brief  description 
of  the  catalyst  layer  models  is  also  presented  with  some  technological  developments  at  nano-scale  to 
improve  the  physio-  and  electro-chemical  properties.  A  test  case  for  a  2D  PEM  cathode  is  also  simulated 
for  different  operating  voltages  to  predict  the  water  saturation  effects. 
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1.  Introduction 

Fuel  cells  have  become  a  pivot  of  energy  research  activities  in 
the  present  decade.  With  increasing  energy  demands  and  depleting 
organic  fuels,  a  need  for  sustainable  and  efficient  energy  produc¬ 
tion  had  never  been  felt  so  urged  as  of  today.  With  many  different 
alternative  proposals  provided  by  the  scientific  and  engineering 
communities,  fuel  cells  stand  a  biased  position  because  of  their 
high  efficiency  with  a  byproduct  of  low  to  zero  greenhouse  gas 
emissions  and  abundance  of  fuel  availability.  Among  many  differ¬ 
ent  types  of  fuel  cells,  polymer  electrolyte  membrane  (PEM)  fuel 
cells  have  taken  the  lead  because  of  their  low  operating  parame¬ 
ters,  cost  effectiveness,  high  current  density  and  compactness  for 
mobile  applications  [1-3]. 

An  outlook  of  PEM  fuel  cells  has  a  deceptive  presentation  as 
being  very  simple  and  straightforward  piece  of  equipment  in  both 
making  and  service.  However,  turning  around  the  coin  indicates 
that  they  are  not  more  simpler  than  any  other  energy  production 
devices,  and  quantifying  and  measuring  all  the  processes  and  phe¬ 
nomena  inside  PEM  fuel  cells  is  not  only  impossible  [3]  but  the 
highly  reactive  environment  also  makes  it  quite  difficult  to  mea¬ 
sure  even  simple  parameters  like  temperature,  pressure,  electric 
potential  and  species  gradients,  etc.  [4].  In  recent  years,  much  criti¬ 
cal  work  has  been  performed  in  various  disciplines  of  PEM  fuel  cells 
from  basic  electro-chemistry  to  design  of  stacks,  but,  numbers  of 
issues  are  still  pending  and  need  to  be  resolved  for  commercial  via¬ 
bility  and  many  improvements  are  deemed  necessary  to  remove 
the  big  question  mark  about  the  future  of  PEM  fuel  cells  as  an 
alternative  energy  production  unit. 

It  has  been  well  established  that  cathode  performance  is  one 
of  the  key  issue  still  under  intensive  investigation  without  any 
proper  remedy  yet  proposed  [5].  The  important  factors  affecting 
the  cathode  performance  are  [6]; 

•  slow  reaction  kinetics, 

•  formation  of  liquid  water  and  water  management, 

•  thermal  management. 

In  PEM  fuel  cells,  the  oxygen  reduction  reaction  (ORR)  is  the  rate 
determining  step  for  the  overall  electro-chemical  reaction.  Despite 
the  active  research  in  improving  the  physio-chemical  behavior 
of  the  cathode  catalyst  it  has  been  determined  that  the  ORR  is 
about  four  to  six  times  slower  to  the  hydrogen  oxidation  reactions 
(HOR)  occurring  at  the  anode  [6,7].  Formation  of  liquid  water  at 
the  cathode  of  PEM  fuel  cell  is  an  another  major  contributor  to  the 
under-grade  performance  of  the  cathode  especially  at  high  loads 
by  blanketing  the  reaction  sites  by  making  them  unavailable  for 
three-phase  contact. 

Different  processes  contributing  to  water  formation  or  removal 
during  the  operation  of  PEM  fuel  cells  at  the  cathode  are;  [8]  (the 
negative  mechanisms  in  water  source  represent  removal  of  water 
content  while  increase  in  water  quantity  inside  the  fuel  cell  is  rep¬ 
resented  by  the  positive  sources); 


Nomenclature 

Aagg  Effective  agglomerate  surface  area  (m2  nr3) 

C^lk  Local  02  concentration  (mol  nrr3 ) 

Db,  Dj  Bulk  diffusivity  of  species  i  (m2  s-1 ) 

Dk  Knudsen  diffusivity  (m2  s-1 ) 

Deff  Effective  diffusivity  of  species  (m2  s-1 ) 

F  Faraday’s  constant 

hv  Interstitial  heat  transfer  coefficient  (W  nrr3  K-1 ) 

K  Permeability  (m2) 

k  Reaction  rate  constant  (s-1 ) 

/<con  Condensation  rate  constant  (s-1 ) 

/<ev a  Evaporation  rate  constant  (Pa  s-1 ) 

kn  Knudsen  number 

l  Characteristic  flow  dimension  (m) 

i  Current  density  (A  nrr 2 ) 

Mi  Molecular  weight  of  species  (kg  mol-1 ) 

Pc  Capillary  pressure  (Pascals) 

R  Universal  gas  constant  (J  mol-1 1<-1 ) 

r  Radius  (m) 

X  Species  mass  fraction 

z  Number  of  electrons  consumed  per  mole  of  reactant 

Greek  Letters 

(Sfiim  Thickness  of  electrolyte  film  covering  an  agglomer¬ 
ate  (m) 

£agg  Proportion  of  electrolyte  in  agglomerate 
s  Porosity  of  material 

0  Theile’s  modulus 

X  Mean  free  path  of  molecules  (m) 

p  Density  (kg nrr3) 

a  Surface  tension  (N  nrr1 ) 

Subscripts  and  superscripts 
agg  Agglomerate 

c  Catalyst  layer 

eff  Effective 

/  Fluid  phase 

i  Species 

Pt  Platinum 

s  Solid  phase 

v  Void  space 


•  oxygen  reduction  reaction  (positive) 

•  electro-osmotic  diffusion  (positive) 

•  condensation  of  water  vapors  (positive) 

•  back  diffusion  (negative). 

•  evaporation  (negative). 

The  ion  transport  in  form  of  H30+  uses  water  molecules  as  a  car¬ 
rier  from  anode  to  cathode  of  a  PEM  fuel  cell.  It  is  estimated  that 
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Fig.  1.  Comparison  of  current  density  at  different  voltages  with  various  levels  of 
water  flooding  in  the  porous  media. 


one  to  five  molecules  of  water  are  dragged  per  proton  migration 
from  anode  to  cathode  side  [9,10].  Similarly,  along  the  production 
of  water  due  to  ORR,  condensation  of  water  vapors  also  proves  to 
be  handy  in  formation  of  liquid  water  if  the  vapor  quantity  exceeds 
its  saturation  limit  in  the  inlet  air  supply  to  PEM  fuel  cell.  The  anti¬ 
measures  for  the  formation  of  water  are  back  diffusion  resulting  due 
to  the  concentration  difference  in  water  across  the  anode  and  the 
cathode,  and  the  evaporation  of  liquid  water  due  to  high  inlet  tem¬ 
perature  or  saturation  pressure.  If  the  formation  rate  of  water  is  not 
balanced  by  the  removal  rate,  accumulation  of  liquid  water  occurs 
at  the  cathode  resulting  in  water  flooding.  This  non-equilibrium  of 
production  and  removal  is  known  to  cause  major  performance  hold 
ups  to  PEM  fuel  cells  in  terms  of  efficiency,  stability  and  reliability 
[11,12]. 

Although,  water  formation  has  been  labeled  as  one  of  the  perfor¬ 
mance  defectors  in  PEM  technology,  many  processes  inside  the  PEM 
fuel  cells  are  itself  highly  water  dependent.  As  already  stated,  the 
proton  migration  from  anode  to  cathode  i.e.,  the  protonic  conduc¬ 
tivity  of  the  membrane  material  incorporated  in  low  temperature 
(<100  °C)  PEM  fuel  cells,  is  highly  water  dependent.  The  dryness  of 
the  membrane  will  render  it  from  low  to  zero  conductivity  caus¬ 
ing  major  suffering  in  performance  by  considerably  increasing  the 
ohmic  losses  [10].  To  ensure  proper  hydration  of  the  membrane, 
a  balance  between  inlet  humidification  and  evaporation  rate  has 
to  be  maintained.  So,  overall,  the  formation  and  removal  rate  of 
water  has  to  be  closely  monitored  and  balanced  not  only  to  avoid 
flooding  of  the  cathode  but  also  for  proper  wetting  of  the  mem¬ 
brane. 

As  discussed  above,  the  water  management  problem  is  one  of 
the  major  issues  related  to  the  optimum  performance  and  sta¬ 
ble  operation  in  PEM  fuel  cells.  A  comparison  of  different  water 
flooding  levels  is  given  in  Fig.  1  for  the  current  produced  for  each 
operating  voltage.  As  it  can  be  seen  that,  when  the  water  flood¬ 
ing  increases,  the  maximum  current  density  for  a  specific  voltage 
decreases  considerably.  At  low  current  densities,  the  current  den¬ 
sity  is  almost  the  same  for  all  levels  because  at  such  low  operating 
conditions  the  reaction  rate  is  quite  low  and  water  formation  due  to 
both  ORR  and  electro-osmatic  drag  is  not  significant  but  at  higher 
levels  of  water  flooding  there  is  decrease  of  almost  80%  in  the  cur¬ 
rent  density  produced  at  0.4  V. 

Many  reviews  of  PEM  fuel  cell  modeling  have  been  published 
by,  e.g.,  Biyikoglu  [13],  Cheddie  and  Munroe  [4],  Wang  [14],  Har- 
aldsson  and  Wipke  [15],  Siegel  [3]  and  Mench  [16]  etc.  Most  of  the 
reviews  were  conducted  for  general  models  related  to  conserva¬ 
tion  equations,  spatial  dimensions  and  level  of  model  complexity. 
The  present  work  is  limited  to  two-phase  flow  models  as  the  water 
management  still  remains  one  of  the  key  issues  for  PEM  fuel  cells. 
Also,  a  brief  insight  will  be  provided  for  micro-scale  model  devel¬ 


opments  in  PEM  fuel  cell  both  in  terms  of  Computational  Fluid 
Dynamics  (CFD)  analysis  and  catalyst  layer  modeling. 

2.  Classification  of  models 

Because  of  vast  diversity  of  technologies  incorporated  in  a  sin¬ 
gle  PEM  fuel  cell,  it  is  quite  difficult  to  classify  and  fix  a  certain 
model  to  a  particular  subclass,  e.g.,  even  in  CFD  modeling,  equa¬ 
tions  of  voltages  need  to  be  solved  for  electro-chemical  reaction 
rates  besides  verification  of  modeling.  Similarly,  electrical,  heat  and 
species  transport  losses  need  to  be  accounted  for  the  simplest  of  the 
models.  Apart  from  this,  anisotropy  of  material  properties  extend 
the  models  floating  in  various  domains. 

Many  authors  have  already  attempted  to  classify  PEM  fuel  cell 
models  according  to  their  own  dominions,  e.g.  Khan  [8]  has  clas¬ 
sified  PEM  models  based  on  thermal  analysis  (isothermal  and 
non-isothermal),  flow  (single-  or  two-phase)  and  the  electro¬ 
chemical  model  used  to  simulate  the  reactions  in  the  catalyst  layer. 
Similarly,  Cheddie  and  Munroe  [4]  have  categorized  PEM  fuel  cell 
models  based  on  modeling  approach  used,  i.e.,  analytical  models, 
semi-empirical  models  and  mechanistic  models.  Analytical  mod¬ 
els  represent  the  simplest  of  all  as  many  simplifying  assumptions 
are  employed  that  results  in  approximate  analytical  voltages  versus 
the  current  density  relationships  [17-19].  In  semi-empirical  mod¬ 
els,  empirically  determined  properties  are  used  with  theoretical 
differential  and  algebraic  equations  [20],  while,  mechanistic  mod¬ 
els  have  been  more  popular  in  modeling  in  which  the  differential 
equations  are  derived  from  physics  and  electro-chemistry  of  the 
internal  phenomena  in  PEM  fuel  cells  [21-23]. 

Siegel  [3]  has  categorized  PEM  fuel  cell  models  in  his  review  arti¬ 
cle  based  on  geometric  constraints  of  the  models  from  one  to  three 
dimensions.  Classification  based  on  the  length  scales  of  the  compu¬ 
tational  domains  has  also  been  proposed  by  Mench  [16]  and  Djilali 
and  Sui  [24].  The  length  scale  varies  from  the  molecular  level  to  full 
system  size  with  different  purposes  and  outcomes.  The  molecular 
models  deal  more  with  an  attempt  to  model  transport  of  charge, 
mass  and  heat  to  interpret  the  limitations  that  significantly  affect 
the  overall  performance  of  fuel  cells  [16,25-28]  whereas,  system 
or  stack  models  deal  more  with  efficiencies,  losses  and  geometric 
limitations  of  the  complete  energy  system  [19,29-38]. 

3.  Macroscopic  models  and  challenging  issues 

Basically,  a  fuel  cell  is  an  electro-chemical  device  that  converts 
the  chemical  energy  into  the  electrical  energy  without  any  inter¬ 
mediating  assistance  or  device.  Main  components  of  a  single  PEM 
fuel  cell  can  be  listed  as; 

•  membrane 

•  catalyst  layer 

•  gas  diffusion  layer  (GDL,  also  known  as  porous  transport  layer, 
PTL) 

•  current  collectors 

•  flow  channels. 

Membrane  is  the  parting  component  between  anode  and  cath¬ 
ode  sides  of  the  fuel  cell  while  catalyst  layer,  GDL  and  current 
collectors  on  either  side  of  the  membrane  constitute  the  electric 
poles  of  the  cell.  Fuel  (hydrogen)  is  fed  to  the  anode  side  of  the  fuel 
cells,  and  is  distributed  on  the  catalyst  layer  by  the  GDL  to  produce 
electrons  and  protons  according  to  Eq.  (1); 

2H2  ->  4H+  +  4e_  (1) 

Electrons  are  forced  to  flow  through  the  external  path,  while, 
protons  migrate  internally  through  membrane  by  selecting  such  a 
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Table  1 

Governing  equations  for  PEM  fuel  cell  models  and  applicable  component  of  PEM  fuel  cell. 


Equation 

Region  of  application 

Remarks 

1. 

Continuity 

V  ■  (pv)  =  Sm 

CL,  GDL,  flow  channels 

X 

2. 

Momentum 

V  ■  (pvv)  = 

-Vp  +  V  •  (?)  +  pg  +  F 

CL,  GDL,  flow  channels 

X 

3. 

Species  transport 

v.  (p?Xf)  =  -V  ■  j(  +  J?i  +s( 

CL,  GDL,  flow  channels 

X 

4. 

Energy  equation 

V  •  ( pcpJ} )  =  V  •  (keffVTf)  + 

Sf0  =  V  •  (/ce,fV7})  +Ss 

All  regions 

One  equation  or  LTNE  approach 

5. 

Electric  potential 

-V  ■  (orsV0s)  = 

S(/)s—  V  ■  (CTmV0m)  —  S</>m 

All  regions  except  flow  channels 

Solid  phase  and  membrane  phase 
potential 

6. 

Secondary  phase 

V  ■  (p/Vs)  =  Rw 

All  regions  except  current  collectors 

Multi-phase  models  in  gas 
channels,  water  saturation 
equation  in  CL  and  GDL,  water  flux 
in  membrane 

material  that  poses  high  electron  resistivity  and  proton  conductiv¬ 
ity  (the  detailed  structure  of  membrane  materials  can  be  found  in 
[39]).  Oxygen  and  charged  entities  (e_  and  H+)  transported  from 
anode  combine  at  the  cathode  to  produce  water  as  a  product. 


•  species  transport 

•  heat  conduction 

•  electrical  conduction 

•  water  saturation. 


02+4H++4e-  -*  2H20  (2) 

Eqs.  (1)  and  (2)  represent  the  half-cell  reactions  of  anode  and 
cathode  sides  respectively  and  are  catalyzed  with  platinum  present 
in  the  catalyst  layer  on  either  sides. 

With  advances  in  computer  technologies  and  enhanced  speeds, 
CFD  modeling  approach  has  provided  scientists  and  engineers  with 
internal  anatomy  of  fundamental  processes  of  a  PEM  fuel  cell  [8].  To 
date  many  models  for  PEM  fuel  cells  have  been  proposed  with  vary¬ 
ing  complexity  and  texture,  however,  there  is  no  single  complete 
model  that  would  effectively  and  efficiently  explain  and  chalk  out 
the  phenomena  altogether.  Nevertheless,  for  a  descriptive  model  of 
PEM  fuel  cell,  Biyikoglu  [13]  has  outlined  basic  conditions  or  pro¬ 
cesses,  inclusion  of  which  can  result  in  a  much  descriptive  model, 
given  as; 

•  balanced  current  distribution 

•  control  of  water  flow 

•  efficient  removal  of  liquid  water 

•  removal  of  excessive  heat. 

Although,  basic  outline,  as  explained  above,  is  very  useful  in 
developing  a  CFD  model  of  PEM  fuel  cells,  still  ultimate  and  com¬ 
plete  PEM  fuel  cell  model  is  quite  difficult  to  achieve  due  to  inherent 
limitations  of  the  analysis  and  outputs  desired.  The  main  limita¬ 
tions  still  blocking  the  researchers  from  attaining  the  ultimate  goal 
of  completeness,  as  given  by  Mench  [16],  are  the  inclusion  of  the 
physico-chemical  phenomena,  knowledge  of  the  transport  phe¬ 
nomena,  computational  power  and  proper  validation  of  the  models. 

The  governing  equations  used  in  PEM  fuel  cell  models  are  given 
in  Table  1  with  applicable  component  regime  and  the  source  terms 
for  each  governing  equations  are  given  in  Table  2  for  a  typical  PEM 
fuel  cell  CFD  model.  The  detailed  description  of  the  governing  equa¬ 
tions  along  with  nomenclature  and  source  term  evaluation  for  all 
the  fundamental  processes  can  be  found  in  [8],  here  some  of  the 
important  and  usually  ignored  factors  are  outlined. 

3.1.  Anisotropy  of  physical  properties 

Regarding  the  porous  media  of  PEM  fuel  cells,  it  is  revealed  that 
it  comprises  of  fibrous  media  that  has  significant  anisotropy  due  to 
its  orientation  of  the  fibers.  Due  to  this,  the  in-plane  and  through- 
plane  properties  vary  significantly  [40,41].  The  major  properties 
influenced  by  anisotropy  are; 


For  species  transport  in  the  porous  media,  most  of  the 
authors  use  Bruggeman’s  correction  factor.  However,  accounting 
for  anisotropy,  the  effective  diffusion  coefficient  is  a  tensor  rather 
than  scalar  quantity.  As  presented  by  Nam  and  Kaviany  [42],  the 
effective  diffusion  coefficient  in  the  porous  media  in  PEM  fuel  cells 
is  better  depicted  by  using  percolation  theory,  given  as; 


Df=/(e)  X  D‘g 

*•>-•(£?)’ 


0.521  in-plane 
0.785  through-plane 


(3) 


where  ev  is  the  percolation  critical  value  and  has  been  reported 
to  be  0.11  and  0.13  by  Pharoah  et  al.  [43]  and  Liu  and  Wang  [44], 
respectively.  The  results  produced  by  the  anisotropic  diffusion  coef¬ 
ficient  reveal  that  the  gas  flow  is  much  higher  in  in-plane  direction 
than  through-plane  directions  suffocating  the  reaction  sites  in  the 
catalyst  layer.  Similarly,  the  thermal  conductivity  of  the  porous 
media  is  almost  14  times  more  in  in-plane  direction  and  tempera¬ 
ture  differential  of  more  than  5  °C  in  through-plane  was  observed 
in  work  by  Pasaogullari  et  al.  [41].  Meng  [40]  has  also  suggested 
anisotropic  electrical  conductivity  of  the  porous  media  and  a  dif¬ 
ference  of  more  than  2000  (s  m-1 )  was  reported  to  exist  in  in-plane 
and  through-plane  directions.  Similarly,  apart  from  anisotropic 
effects,  the  experimental  and  CFD  modeling  results  also  vary  as  the 
temperature  gradients  observed  are  much  lower  while  simulated 
numerically  using  the  parallel  resistance  approach,  given  as  [43]; 


keff  =  skv  +  (!  -£)1<s 


(4) 


Table  2 

Source  terms  for  governing  equations. 


Equations 

Source  term 

Gas  diffusion  layer  (GDL) 

Catalyst  layer  (CL) 

Momentum 

hlORR,02 

0 

hlORR,H20 

0 

^v.i 

f^Phase 

Eq.  (12) 

Eq.  (12) 

i2 

i2  f 

Energy 

ls 

°s,eff 

s  m 

°s,eff  °m,eff 

*7phase 

^phase  x  tlfg 

th  phase  x  tlfg 

< lit 

hv(Ts  -  Tf) 

0 

Qorr 

0 

(0m-0s)x  V  i 

Charge 

S<Ps 

0 

Vi 

S<pm 

0 

-Vi 
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where  ks  and  kv  are  solid  and  void  thermal  conductivities,  respec¬ 
tively.  The  conduction  of  heat,  as  predicted  by  Eq.  (5)  is  assumed 
to  travel  either  through  void  or  solid  regions  contrary  to  combined 
path.  Pharoah  et  al.  [43]  suggested  a  novel  approach  by  stacking 
the  solid  and  void  regions  so  that  all  of  the  heat  is  passed  through 
combined  solid  and  void  regions  termed  as  the  network  resistance 
approach. 

keff  =  eks  +  (1-  Wv  (5) 

The  heat  transport  predicted  by  utilizing  the  thermal  coefficient 
as  given  by  Eq.  (5)  is  approximately  the  same  as  encountered  in 
reality  [43]. 

Similarly,  most  of  the  two-phase  models  use  Leverett  function 
to  simulate  the  water  flow  in  PEM  fuel  cells  [8,41 ,45].  This  function 
is  highly  dependent  on  the  permeability  of  the  porous  media  which 
in  turn  is  actually  anisotropic  in  nature,  as  given  by  [41  ]; 

jx  =  Qkxx  and  jy  =  C2  kyy  (6) 

where  jx  and  jy  are  the  mass  flux  of  liquid  water  in  respective 
directions,  kxx  and  kyy  represent  the  anisotropic  permeability  in 
the  principle  directions  of  the  porous  media.  It  is  found  that  water 
saturation  levels  are  always  higher  in  this  case,  thus  reducing  the 
overall  performance  of  the  PEM  fuel  cells  [41  ]. 


3.2.  Local  thermal  non- equilibrium  (LTNE)  approach 


The  condition  to  use  one  equation  model  for  heat  transport  or 
local  thermal  equilibrium  (LTE)  is  only  valid  when  the  tempera¬ 
ture  difference  between  the  fluid  and  the  solid  phases  of  the  porous 
electrode  is  much  lower  than  the  overall  system  temperature  differ¬ 
ence.  Since,  PEM  fuel  cells  are  low  operating  temperature  devices, 
the  magnitude  of  both  the  temperature  differences  between  phases 
and  the  overall  temperature  difference  is  almost  same,  given  math¬ 
ematically  as  [46]; 


O 


^^loc 

ATS  yS 


(7) 


where  zATloc  is  local  temperature  difference  between  the  phases 
while  ATsys  represents  the  overall  system  temperature  differ¬ 
ence.  In  LTNE  approach  both  the  phases  i.e.,  void  and  solid 
portions  of  the  porous  media  are  represented  by  separate  equa¬ 
tions  and  interlinked  through  a  volumetric  heat  transfer  coefficient 
hv  (W  m-3 1<-1 ).  The  exact  value  of  hv  has  not  been  measured  yet 
but  experimental  results  obtained  for  the  aluminum  foam  rang¬ 
ing  from  3  x  104  to  1.5  x  105  (Wnrr3  K_1)  [46]  has  been  used  by 
some  authors.  To  date  many  researchers  have  implemented  LTNE 
approach  in  both  single  phase  and  two-phase  flow  regimes  but  lim¬ 
iting  the  geometry  to  2D  only  [8,23,46-49].  Since  many  parameters 
in  PEM  fuel  cell  are  temperature  dependent,  LTNE  approach  to  3D 
models  need  to  be  evaluated  and  compared  with  the  local  thermal 
equilibrium  or  the  so-called  one  equation  models. 


3.3.  Knudsen  diffusion 

Usually  effective  diffusion  coefficient  modified  by  Bruggeman’s 
correction  is  employed  in  species  transport  of  oxidants  and  fuels 
in  the  porous  media  of  the  PEM  fuel  cells.  However,  the  texture 
of  porous  media  is  very  complex  and  the  relative  influence  of 
ordinary  diffusion  or  Knudsen  diffusion  on  species  transport  is  gov¬ 
erned  by  the  pore  geometry  [50].  According  to  Malek  and  Coppens 
[51],  for  the  media  with  pore  dimensions  of  2-50  nm,  Knudsen 
diffusion  is  the  predominant  transport  mechanism  which  results 
from  the  collision  of  gas  molecules  with  the  pore  walls  instead 


of  intra-molecular  collision  (Brownian  motion).  The  Knudsen  dif¬ 
fusion  coefficient  for  CFD  calculations  can  be  used  in  the  form  of 
[52]; 


The  effective  diffusion  coefficient  based  on  both  bulk  diffusion 
and  Knudsen  diffusion  can  be  calculated  as  [53]; 


where /[£)  is  the  correction  factor  for  the  porous  medium  and  can 
be  evaluated  according  to  Eq.  (3). 


3.4.  Modeling  software  and  solutions 

With  in-house  self-developed  and  open-source  CFD  models, 
many  commercial  software  products  are  available  in  the  market 
that  have  proved  to  be  very  efficient  and  robust.  Regardless  of  inher¬ 
ited  disadvantage  of  limited  freedom  in  equation  manipulation  and 
controls,  many  researchers  have  opted  for  commercial  software 
products  as  prime  CFD  tool.  Among  many,  the  most  commonly 
used  are  FLUENT®,  COMSOL®,  STAR-CD®  and  CFD-ACE+®,  and  the 
contribution  of  each  in  CFD  modeling  is  shown  pictorially  in  Fig.  2 
[3]. 

ANSYS®  Fluent®  [54]  is  a  powerful  commercial  software  avail¬ 
able  in  the  market  offering  sophisticated  numerics  and  robust 
formulations  including  a  pressure-based  segregated  and  coupled 
solvers,  and  a  density  based  solver  technique  to  ensure  optimum 
and  reliable  results.  It  is  well  suited  for  a  number  of  complex  physi¬ 
cal  models  utilizing  unstructured  meshes  both  for  2D  and  3D  cases 
based  on  finite  volume  method  (FVM).  Meshes  can  be  created  using 
ANSYS®  meshing  products  or  other  robust  products  like  ICEM®  and 
GAMBIT®.  To  further  enhance  the  flexibility  of  model  variant  situ¬ 
ations,  Fluent®  provides  the  use  of  user-defined-functions  (UDFs) 
that  help  to  tailor  the  model  to  specific  needs  or  requirements.  With 
the  increase  in  demands  for  fuel  cells  CFD  analysis,  an  add-on  is 
provided  to  simulate  typical  PEM  fuel  cells.  In  this  module,  the  two- 
phase  flow  in  the  porous  media  is  solved  using  the  Leverett  function 
while  in  gas  channels  it  is  assumed  that  both  gases  and  water  flow 
with  the  same  velocity  in  form  of  fine  mist.  The  potential  equations 
for  both  solid  and  membrane  phases  are  also  solved  to  get  the  inside 
distribution  of  electric  and  ionic  currents,  respectively.  The  addi¬ 
tional  capability  of  Fluent®  software  to  solve  user-defined-scalars 
provides  a  handy  tool  to  simulate  most  of  the  internal  phenom¬ 
ena  occurring  inside  PEM  fuel  cells.  Fluent  also  offers  extended 
post-processing  tool  and  provides  provision  to  export  the  data  to 
other  post-processing  software  solutions.  In  order  to  simulate  the 
electro-chemical  reactions,  a  choice  is  also  provided  to  select  either 
Tafel  formulation  or  the  advanced  Butler-Volmer  kinetics.  Further¬ 
more,  to  enhance  the  flexibility  of  the  PEM  fuel  cell  module,  it  is  also 
made  possible  to  specify  requirement  based  specific  functions  and 
models  through  a  modifiable  source  code. 

COMSOL  Multiphysics®  [55]  is  another  commercial  software 
product  available  widely  for  simulating  fluid  flow  incorporating 
finite  element  method  (FEM)  with  stiff  chemistry  solvers.  This 
software  also  offers  flexible  equation  based  modeling  with  vari¬ 
ous  degrees  of  mesh  complexity.  The  user  can  also  include  partial 
differential  equations  (PDEs)  using  various  formulations  and  the 
software  can  automatically  detect  the  best  possible  solver  and  set¬ 
tings  for  a  particular  problem  along  with  manual  tuning.  Since  fuel 
cells  have  become  a  beacon  of  future  power,  an  add-on  module 
has  been  released  in  the  latest  version  V4.  Compared  to  Fluent® 
in  which  add-on  module  is  only  available  for  PEM  and  solid  oxide 
fuel  cells  (SOFCs),  COMSOL®,  on  the  other  hand,  also  supports  alka¬ 
line,  molten  carbonate  and  direct  methanol  fuel  cells.  As  far  the 
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Fig.  2.  Comparison  of  commercial  and  non-commercial  software  used  for  PEM  fuel 
used  for  3D  modeling.  (For  commercial  software,  data  is  adopted  from  Siegel  [3]). 

mesh  generation  for  COMSOL  applications  is  concerned,  a  built-in 
geometry  and  mesh  generator  is  available  in  standard  package  with 
sound  post-processing  facility. 

STAR-CD®  [56]  and  CFD-ACE+®  [57]  are  other  finite  volume 
based  CFD  software  products  also  commercially  available.  Both 
products  have  strong  parallel  computing  capabilities  with  aptitude 
to  handle  complex  geometries  including  multiphysics  modules  for 
PEM  fuel  cells  that  can  be  handled  quite  easily.  CFD-ACE+®  also 
offers  the  facility  to  carryout  stress  and  strain  analysis  within  the 
PEM  fuel  cells.  Alongside  licensed  commercial  software,  Open- 
FOAM  is  open  code  software  with  pre-developed  modules  for 
simulating  complex  multiphysics  problems  with  good  pre-  and 
post-processing  utilities  also  based  on  finite  volume  technique. 
Additionally,  meshes  from  other  software  can  also  be  imported  and 
manipulated.  The  advantages  offered  by  OpenFOAM  are  that  users 
can  extend  or  create  own  libraries  and  manipulate  the  solver  to  suit 
the  demands.  Although  no  fuel  cell  modules  are  provided  built-in 
but  with  good  coding  knowledge,  efficient  and  robust  models  can 
be  developed  that  can  handle  a  variety  of  problems  [3]. 

4.  Two-phase  models 

Water  management  inside  PEM  fuel  cell  is  of  paramount  impor¬ 
tance  because  operation  of  a  PEM  fuel  cell  is  highly  dependent  on 
the  water  content  as  protonic  conductivity  of  the  membrane  mate¬ 
rial,  typically  used  for  low  operating  temperature  fuel  cells,  is  highly 
dependent  on  the  amount  of  water  present  [39,58].  Decrease  in 
water  content  can  cause  dry  out  of  membrane  and  reduce  protonic 
migration,  thus,  it  is  very  essential  that  the  membrane  remains  fully 
hydrated  all  the  times.  However,  excessive  water  inside  PEM  fuel 
cell  can  lead  to  clogging  of  channels  [59,60],  and,  flooding  of  the 
catalyst  layer  [59,61,62]  and  long  term  liquid  water  accumulation 
inside  PEM  fuel  cell  is  also  one  of  the  major  contributor  to  the  degra¬ 
dation  of  the  catalyst  and  its  carbon  support  material,  ionomer 
poisoning  and  hydrophobicity  loss  [63].  In  this  section,  an  insight 
into  selected  two-phase  models  is  presented  with  a  summarized 
outline  in  Table  3.  The  list  presented  here  is  not  an  exhausted  one 
but  it  has  been  made  to  include  different  variety  of  two-phase  CFD 
models  used  in  the  research  society  in  last  10  years. 

4.1.  CFD  modeling  review 

He  et  al.  [64]  performed  multi-phase  simulations  of  PEM  fuel 
cell.  In  their  model,  the  droplet  size  of  water  was  taken  as  the  prime 
parameter  to  study  the  multi-phase  flow.  A  multi-phase  mixture 
model  was  applied  to  a  3D  geometry.  The  interaction  effect  between 
the  phases  was  studied  by  considering  droplet  size,  drag  coefficient, 
Reynolds  number,  velocity  and  droplet  relaxation  time.  The  equa¬ 
tion  for  calculating  the  droplet  size  was  adapted  from  Zhang  et  al. 
[60].  The  advantage  associated  with  this  model  is  that  it  includes 
the  effect  of  liquid  water  removal  from  the  channels  of  PEM  fuel  cell. 


CFD.  (a)  Trend  of  software  used  for  CFD  modeling  in  overall;  (b)  trend  of  software 


In  this  work  all  the  simulations  were  performed  on  a  commercial 
CFD  software. 

Yuan  and  Sunden  [65]  also  developed  a  3D  model  to  under¬ 
stand  the  effect  of  liquid  water  on  cell  performance.  The  model 
presented  is  a  half-cell  model  considering  cathode  only  because  of 
its  slower  kinetics.  The  electro-chemical  reactions  were  modeled 
to  occur  in  a  thin  layer  while  simulating  the  flow  in  GDL  and  the 
channels.  The  salient  features  of  this  work  were  the  use  of  com¬ 
bined  thermal  boundary  conditions  and  mass  transfer  along  with 
the  effect  of  saturation  on  current  density  profile.  The  calculation 
domain  was  discretized  by  finite  volume  method  and  a  combina¬ 
tion  of  uniform  and  non-uniform  grid  spaces.  The  simulations  were 
performed  using  an  in-house  CFD  code.  The  saturation  was  evalu¬ 
ated  based  on  the  saturation  pressure  and  the  local  temperature  of 
the  flow.  Wang  et  al.  [66]  also  proposed  a  two-phase  model  of  PEM 
fuel  channels  to  simulate  the  flow  of  liquid  water  and  gaseous  reac¬ 
tants.  Darcy’s  law  and  multi-phase  mixture  model  was  employed 
for  estimation  of  the  key  parameters  such  as  overall  pressure  drop 
and  liquid  saturation  profile  along  the  channel  and  flow  analogy 
to  random  porous  media.  The  physical  model  used  in  this  work 
consisted  of  a  single  straight  channel  for  the  cathode  side.  Fully 
or  partially  humidified  air  feed  was  used  at  the  inlet  and  the  water 
produced  as  a  result  of  electro-chemical  reactions  was  injected  into 
the  channel  along  its  length. 

He  et  al.  [21]  performed  a  2D  analysis  of  a  PEM  fuel  cell  with 
interdigitated  flow  field  by  solving  two-phase,  multi-component 
transport  equations.  Darcy  law  was  used  to  simulate  the  gas  phase 
transport,  while  for  liquid  water,  both  the  capillary  pressure  and  the 
shear  force  between  the  two-phase  was  considered  as  the  prime 
transport  mechanisms.  The  modeled  region  consisted  of  a  GDL  and 
the  current  collector.  The  electro-chemical  reactions  were  assumed 
to  occur  at  the  boundary  of  GDL,  i.e.,  a  thin  layer  model  was  used. 

Le  and  Zhou  [67]  presented  a  general  model  of  PEM  fuel  cell 
which  accounted  for  fluid  flow,  heat  transfer,  species  transport, 
electro-chemical  reactions  and  water  saturation.  Detailed  thermo¬ 
electro-chemistry  studies  were  carried  out  on  a  3D  geometry  where 
saturation  effects  were  evaluated  with  explicit  gas-liquid  inter¬ 
face  tracking  using  VOF  multi-phase  model.  Commercial  software 
Fluent®  was  employed  for  simulations  and  calculations.  In  this 
model,  all  the  components  of  a  complete  single  fuel  cell  unit  were 
included  to  broaden  the  results  scope  and  the  effect  of  liquid  sat¬ 
uration  along  with  porous  media  was  also  studied.  The  flow  field 
design  employed  in  this  model  was  the  serpentine  flow  field.  Ini¬ 
tially,  water  droplets  of  0.4  mm  were  assumed  suspended  in  the 
flow  field  and  their  behavior  was  studied  at  different  time  steps  at 
an  operating  voltage  of  0.5  V. 

Hwang  [23]  has  also  developed  a  2D  model  of  PEM  fuel  cell 
in  which  two  separate  momentum  equations  were  applied  for 
gaseous  and  liquid  flow.  Heat  distribution  was  simulated  using  the 
LTNE  approach  that  considered  separate  energy  equations  for  fluid 
and  solid  components  in  the  GDL.  Also,  irreversible  heat  generation 
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Table  3 

Summarized  features  of  the  review  articles. 
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Author 


1.  Heetal.  [64] 


2.  Yuan  and  Sunden 
[65] 


3.  Heetal.  [21] 


4.  Le  and  Zhou  [67] 


5.  Hwang  [23] 


Chemical  Saturation  equation  Model  Geometry  Major  conclusions 

kinetics  validation 


Butler-Volmer  |  ( spis )  +  V  •  (p, ^  |  (Pc)  Vs)  =  S 


Tafel  equation  s  =  p^w  pns<j!wv 


Tafel  Equation  vl  =fvg  -  DcVs 


Butler-Volmer  ^  (£S/p,)  +  V  (stpiv)  =  S 


Butler-Volmer 


\  V  ■  ( pugug )  =  -e  [  Vpc  +  (pw  -  pg)g] 

,  V-7  /  T — 7  S  SflWUW 

+v  ■  (/rwVuw) - ; - 


Experimental  3D 
Fluent  PEM 
module 


X  3D 

Cathode  only 


Experimental  2D 

Cathode  only 


Visual  3D 


Large  water  droplets 
increase  saturation 
Liquid  water  hinders 
heat  transfer  in  gas 
diffusion  layer  and 
catalyst  layer  causing 
hot  spots 
Co-flow  pattern  is 
disadvantageous  for 
PEM  fuel  cells  as 
compared  to 
counter-flow  pattern 
The  non-uniformity  in 
the  current  distribution 
resulted  because  of 
uneven  distribution 
and  transport  of  the 
reacting  air 
Higher  saturation 
levels  resulted  in  low 
current  densities 
Low-operating 
parameters 
(temperature  and 
humidity)  also 
decreases  the  current 
density 

In  interdigitated  flow 
fields  along  with 
evaporation  of  liquid 
water,  the  liquid  water 
transport  also  acts  as 
water  removal 
mechanism 
Higher  oxygen  flow 
improves  the  reactant 
transport 
The  thickness  of 
electrodes  effects  the 
current  density 
For  interdigitated  flow 
fields,  a  shorter  rib 
dimensions  improves 
the  performance 
Liquid  droplets  cause 
high  pressure  drop 
regions 

The  turning  areas  of 
the  serpentine  flow 
field  showed  increased 
water  droplet 
formation  causing  a 
substantial  air  flow 
blockage 

The  channel  structure 
significantly  influences 
the  water  distribution 
inside  the  cell 
High  saturation  causes 
excessive 

concentration  losses 
Area  under  the 
channels  give  higher 
reaction  rates 
However,  the  presence 
of  liquid  water  reduced 
the  overall 

temperature  of  the  cell 


Experimental  2D  Increasing  the  current 

collector  temperature 
reduces  water 
saturation  by  helping 
evaporation. 
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Table  3  ( Continued ) 

Model  Geometry  Major  conclusions 

validation 


Author  Chemical  Saturation  equation 

kinetics 


6.  Senn  and 

Poulikakos  [45] 


Tafel  equation 


cos  6 


■J(s) 


Pc  =  ■ 

J(s )  =  1.417s  -  2.12s2  +  1.263s3 


X 


2D 


7.  Siegel  et  al.  [73]  Agglomerate  u  ■  Vs  =  V  ■  (d^  Vs)  +  S 


Experimental  2D 


8.  Yu  et  al.  [74]  Butler-Volmer  V  ■  ( epuCa )  =  -V  Ja+S 


Experimental  3D 


9.  You  and  Liu 
[69,70] 


Butler-Volmer 


(pC“)  +  V  ■  ( yapuCa )  =  V  •  (spDVC“)  + 
at 


V- 


*E  [PkskD«  (vq 


VC“)] 


-  V- 


Experimental  2D 


Decreasing  the 
dimensions  of  current 
collector  width, 
channel  width  and  GDL 
thickness  improved  the 
performance  of  the 
PEM  fuel  cell 
For  small  channel 
widths,  a  thin  GDL 
performed  better  that 
thick  GDL 
Liquid  water  has 
significant  effect  on 
fuel  cell  performance 
Water  content  in 
catalyst  layer  and 
transport  through 
membrane  also  effect 
the  ohmic  losses  and 
reactant  transport 
properties  in  the 
cathode  catalyst 
20-40%  of  the  liquid 
water  is  accumulated 
in  the  cathode  catalyst 
layer 

Pure  oxygen  instead  of 
inlet  air  increases  fuel 
cell  performance 
A  thin  PEM  fuel  cell  has 
better  performance 
characteristics  than  a 
thicker  one 
High  inlet  velocity 
increases  the  current 
density 

Increasing  current 
collector  width  has  no 
prominent  effect  on 
the  current  density 


The  water  content 
increasing  along  the 
thickness  of  the 
membrane  towards  the 
cathode  indicating 
large  water  migration 
from  anode  to  cathode 
Increasing  the  current 
increases  oxygen 
consumption 
At  the  cathode  side,  the 
protonic  current  is  zero 
at  the  GDL/catalyst 
layer  interface  and 
increases  towards  the 
catalyst/membrane 
interface  and  is  higher 
for  high  current 
density.  Similarly,  the 
overpotential 
corresponds  directly  to 
the  magnitude  of 
current  density 
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Table  3  ( Continued ) 

Author  Chemical  Saturation  equation 

kinetics 

10.  Meng  [80]  Butler-Volmer  v  ■  [- ^  Vs]  -  V  ■  Vp]  =  S 


11.  Coppoet  al.  [75]  Agglomerate  §  +  V  ■  (u/s)  -  V  ■  (e^qVs)  =  S 


12.  Zhou  et  al.  [83] 


Tafel  equation 


^#,  liquid  M  (  kchd  \ 

dx  ~  \RuT#(x)j 


^#,vapor(x) 

J2iN#’iix) 


(P#M-P#,sat(x)) 


13.  Berning  and  Tafel  equation  Schlogl  equation  for  liquid  water  through  membrane 

Djilali  [86] 


Model  Geometry  Major  conclusions 

validation 


Experimental  2D 
and  model 
comparison 


Experimental  3D 


X  2D 


Experimental  3D 


Water  saturation  level 
is  highest  in  the 
catalyst  layer 
Water  saturation  level 
is  minimum  in 
micro-porous  layer 
Micro-porous  layer 
serves  as  a  barrier  for 
liquid  water  and 
prevents  covering  up  of 
the  interface  between 
gas  flow  channel  and 
the  GDL,  thus, 
increasing  the  species 
transport  in  the  porous 
media 

High  temperature 
increases  the  reaction 
rate 

Higher  temperature 
results  in  higher 
membrane 
conductivity 
Higher  temperature 
increases  the  species 
diffusivity 
Higher  water 
diffusivity  results  have 
been  noticed  by 
increasing  the 
temperature  in  a  highly 
hydrophobic  materials 
Water  advection 
increases  at  the 
interface  of  GDL  and 
gas  channel  at  higher 
levels 

Humidification  of 
anode  and  cathode 
inlet  is  very  important 
Water  content  on  the 
anode  side  dominates 
the  membrane 
performance  by 
delivery  the  protons  to 
the  cathode  side 
Liquid  water  injection 
at  anode  improves  cell 
performance  while 
injection  at  cathode 
hinders  the  effective 
water  removal 
Pressure  loss  in  PEMs  is 
one  of  the  major 
parameters  affecting 
the  overall 
performance  of  cell 
High  inlet  humidity  at 
cathode  decreases  the 
overall  system 
performance  by 
increasing  the 
pumping  power 
Distribution  of  three 
dimensional  flow 
velocities,  species 
concentration,  mass 
transfer  rates,  electric 
current  and 
temperature  were 
illustrated 
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Table  3  ( Continued ) 

Author 

Chemical  Saturation  equation 

kinetics 

Model 

validation 

Geometry 

Major  conclusions 

14.  Hu  etal.  [22,87] 

Butler-Volmer  s  =  ^ 

Vp 

Experimental 

3D 

Interdigitated  flow 

fields  enhance  the 
reaction  rates  by 
providing  better  flow 
through  GDL  as 
compared  to  other  flow 
field  distributions 
Water  saturation  levels 
are  low  with 
interdigitated  flow 
fields 

The  humidification  of 
inlet  species  is  more 
significant  when 
utilizing  interdigitated 
flow  fields  as  compared 
to  conventional  flow 
fields 

15.  Wang  et  al.  [66]  X  s  =  {p^/M^~o)-csat  Experimental  2D,  3D  Liquid  water  buildup  is 

2  faster  near  inlet  region 

under  full 

humidification  inlet 
conditions 

Liquid  saturation  level 
upto  20%  was  observed 
near  inlet 

It  was  observed  that 
water  was  trapped  at 
geometric  changes 
(bend  etc.) 


due  to  the  electro-chemical  reactions,  the  ohmic  losses  and  the  heat 
of  evaporation  and/or  condensation  were  explicitly  considered  in 
this  model.  The  calculation  domain  consisted  of  two  porous  layers 
of  the  cathode  of  PEM  fuel  cell.  The  geometry  considered  in  this 
model  was  interdigitated  flow  field  with  two  regions  of  the  porous 
media  distinguished  according  to  the  position,  i.e.,  the  area  under 
the  inlet  or  the  area  under  the  current  collector.  Finite  element 
technique  was  employed  to  the  solution  domain  considering  the 
liquid  flow  force  due  to  the  imbalance  of  the  liquid  water  pressure 
and  the  gaseous  pressure.  The  correlation  provided  by  Leverett  [68] 
was  employed  to  quantify  the  water  saturation  effects  in  the  porous 
media. 

Senn  and  Poulikakos  [45  ]  have  also  presented  a  2D  model  of  PEM 
fuel  cell  that  accounted  for  multi-component  species  diffusion,  for¬ 
mation  of  liquid  water,  heat  transfer  and  electronic  current.  All  the 
governing  equations  in  this  model  were  non-dimensional  and  FVM 
was  employed  to  solve  the  equations.  The  main  focus  of  the  work 
was  to  study  the  effects  of  varying  the  channel,  current  collector 
and  GDL  dimensions.  The  catalyst  layer  was  assumed  to  be  very  thin 
and  treated  as  a  boundary  condition  for  the  electro-chemical  reac¬ 
tions  and  the  heat  source.  The  distinguished  feature  of  this  work 
was  the  introduction  of  a  performance  variable  based  on  the  aver¬ 
age  current  density  to  measure  the  effects  of  mass  transfer,  water 
saturation  and  heat  transfer. 

A  two-phase  flow  model  has  also  been  presented  by  You  and  Liu 
[69,70]  for  the  cathode  of  PEM  fuel  cell.  The  concept  of  multi-phase 
mixture  model  coupled  for  the  porous  media  and  the  gas  chan¬ 
nel  was  implemented  to  study  the  saturation  effects.  The  model 
presented  by  Wang  and  Chen  [71]  was  extended  to  incorporate 
detailed  effects  of  levels  of  multi-phase  mixtures  instead  of  sep¬ 
arate  phases  (i.e.,  two  fluid  model)  including  non-zero  interfacial 
areas.  The  multi-phase  model  used  in  this  work  is  derived  as  given 
by  Wang  and  Chen  [71]  and  Abriola  and  Pinder  [72].  This  model 
was  limited  to  2D  and  the  cathode  of  a  PEM  fuel  cell  and  explicit 
results  were  obtained  for  the  current  density  affected  by  different 
saturation  levels  along  with  inlet  air  humidification. 


Siegel  et  al.  [73]  also  carried  out  comprehensive  modeling  of  a 
2D  PEM  fuel  cell  including  water  transport  within  the  porous  media 
by  the  capillary  pressure.  The  rate  equation  for  electro-chemical 
reactions  was  adapted  to  the  agglomerate  structure  of  the  catalyst 
layer.  Major  focus  of  this  article  was  to  study  the  effects  of  differ¬ 
ent  factors  affecting  the  PEM  performance,  e.g.,  geometry,  porosity, 
and  polymer  properties,  etc.  The  physical  properties  used  in  this 
model  were  derived  by  direct  measurement  from  a  base  case  fuel 
cell  experimental  model.  The  results  obtained  from  these  simula¬ 
tions  were  used  to  study  the  effect  of  liquid  water  on  reaction  rate 
and  local  decrease  in  the  porosity  of  the  catalyst  layer  and  GDL. 

Yu  et  al.  [74]  presented  a  3D  model  with  interdigitated  flow 
field  for  a  PEM  fuel  cell.  In  this  work  two-phase  flow  and  trans¬ 
port  mechanism  was  developed  to  study  the  performance  of  a  cell 
under  different  operational  parameters.  A  detailed  physical  insight 
was  provided  for  the  velocity,  the  species  concentration,  the  water 
content  and  concentration,  and  the  current  density  distribution. 
All  the  model  equations  were  discretized  using  finite  volume  tech¬ 
nique  and  simulated  using  commercial  CFD  software.  Coppo  et  al. 
[75]  also  developed  a  CFD  model  that  incorporates  all  the  major 
phenomena  in  PEM  fuel  cell,  i.e.,  three-phase  flow  (the  third  phase 
refers  to  the  dissolved  phase),  proton  and  species  transport  etc. 
and  agglomerate  model  was  employed  to  determine  the  reaction 
kinetics.  The  pivot  of  this  work  was  to  evaluate  the  temperature 
dependence  of  all  the  physical  properties  used  in  general  PEM  mod¬ 
els.  Moreover,  a  novel  model  was  also  incorporated  to  describe 
the  liquid  water  from  the  porous  media  surface  by  advection  of 
water  droplets  due  to  the  gas  streams  in  the  gas  channel.  All  the 
simulations  were  performed  at  different  temperature  levels  and 
experimentally  verified  for  the  accuracy  of  this  approach. 

Another  complete  2D  PEM  fuel  cell  model  with  two-phase  flow 
was  presented  by  You  and  Liu  [70].  This  model  was  a  continuation 
of  earlier  models  presented  by  same  authors  [69,76].  A  coupled 
flow,  species,  electrical  potential  and  current  density  was  solved 
in  the  flow  channels,  GDLs,  catalyst  layer  and  the  membrane.  The 
coupling  of  governing  equations  on  both  the  cathode  and  the  anode 
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sides  including  the  membrane  provided  a  much  deeper  insight  of 
the  various  parameters  and  water  content.  To  obtain  a  converged 
solution,  the  authors  first  assumed  the  overpotential  at  the  cata¬ 
lyst/membrane  interface  and  calculated  the  local  current  density 
which  in  turn  was  used  to  simulate  the  net  water  flux  across  the 
boundary  and  the  protonic  current.  Finally  the  oxygen  and  the  cal¬ 
culated  water  flux  were  used  as  the  boundary  condition  for  the 
working  domain.  All  the  coupled  equations  were  solved  iteratively 
and  the  average  current  density  was  measured  by  averaging  the 
local  current  densities  along  the  flow  path.  Acosta  et  al.  [77]  also 
presented  a  2D,  non-isothermal,  two-phase  model  for  PEM  fuel  cells 
with  both  conventional  and  interdigitated  flow  field  designs.  A  con¬ 
tinuum  approach  was  utilized  to  simulate  gas  and  liquid  water  with 
extended  Darcy’s  law.  The  physical  properties  used  in  the  model, 
e.g.,  wettability,  permeability  and  porosity  etc.,  were  determined 
experimentally.  For  the  measurements  of  the  saturation  content  in 
the  porous  media,  a  method  based  on  mercury  intrusion  porosime- 
try  was  used  to  quantify  the  capillary  pressure.  This  model  was 
also  limited  to  the  cathode  side  only  and  the  continuum  approach 
was  used  to  solve  the  coupled  governing  equations.  All  the  gov¬ 
erning  equations  were  solved  iteratively  using  an  in-house  code 
called  MUFTE.UG  which  is  based  on  the  concepts  and  algorithms 
presented  by  Flelmig  [78]  or  Class  [79].  For  simulation  of  water  sat¬ 
uration,  the  Eulerian  approach  was  used  where  discontinuities  in 
pressure  between  the  two  fluid  interface  was  balanced  by  employ¬ 
ing  the  capillary  pressure  effects  that  exist  in  the  porous  media. 

Meng  [80]  presented  a  multi-dimensional  two-phase  model 
including  a  micro-porous  layer  sandwiched  between  the  catalyst 
layer  and  the  GDL.  A  micro-porous  layer  has  been  found  to  reduce 
the  water  saturation  levels  thus  enhancing  the  oxygen  transport  to 
the  reaction  sites  and  the  efficiency  [81  ].  This  model  has  enhanced 
techniques  employed  to  properly  incorporate  the  interfacial  liquid 
water  transport  phenomena  between  the  different  porous  media. 
Furthermore,  the  effects  of  the  current  collector  and  the  gas  flow 
channel  on  the  saturation  of  the  porous  media  have  also  been 
explicitly  studied.  The  results  of  this  model  were  verified  with  high- 
resolution  neutron  imaging  data  [82]  and  other  numerical  data. 
Zhou  et  al.  [83]  has  also  presented  a  multi-phase  and  -component 
2D  model  of  a  PEM  fuel  cell  with  pressure  and  phase  change  effects 
to  further  understand  the  influence  of  the  inlet  humidification  and 
pressure.  One  of  the  major  assumptions  in  this  model  was  that  liq¬ 
uid  water  was  assumed  to  exist  in  form  of  small  droplets  with  no 
volume.  Berning  and  Djilali  [84-86]  also  carried  out  series  of  3D 
work  to  study  the  effects  of  temperature  and  water  management 
on  the  performance  of  fuel  cells.  In  their  work,  a  single  fuel  cell 
was  divided  into  one  main  and  three  sub-domains.  All  the  domains 
were  coupled  through  adjustment  of  appropriate  boundary  con¬ 
ditions.  Similarly,  Hu  et  al.  [22,87]  also  developed  and  analyzed  a 
two-phase  PEM  model  considering  species  transport  in  both  anode 
and  cathode  sides.  A  special  consideration  was  given  to  the  impact 
of  ribs  on  the  species  transport  and  SIMPLE  algorithm  with  fourth 
order  Runge-Kutta  method  was  used  in  their  model. 

4.2.  Experimental  analysis  of  water  transport 

Additionally,  many  authors  have  experimentally  studied  water 
transport  for  specific  components  of  PEM  fuel  cells,  i.e.,  the  catalyst 
layer,  the  membrane  and  the  GDL  that  are  of  paramount  impor¬ 
tance  because  individual  component  has  different  behavior  under 
different  saturation  and  water  content  levels.  A  brief  description  of 
such  methods  is  presented  here  with  briefly  stating  the  findings  of 
the  studies. 

There  are  two  types  of  water  transport  in  the  membrane  of  a 
PEM  fuel  cells,  the  electro-osmotic  drag  and  the  back  diffusion  [88]. 
For  complete  understanding  and  accurate  modeling  of  membrane 
materials  it  is  very  essential  to  accurately  estimate  the  electro- 


osmotic  and  back  diffusion  coefficient  of  the  membrane  material. 
Yan  et  al.  [89]  reported  a  value  of  1. 5-2.6  for  drag  coefficient  for 
Nation™  117  for  different  inlet  humidification  conditions.  Simi¬ 
larly  Ge  et  al.  [90]  determined  that  varying  the  thickness  of  the 
membrane  of  the  PEM  fuel  cell  has  significant  effect  on  the  water 
transport  through  the  membrane  and  hydrophobicity  of  the  mem¬ 
brane  material  also  alters  the  absorption  and  desorption  of  water  at 
the  membrane/catalyst  interface  thus  impacting  the  overall  water 
transport  through  the  membrane  material. 

The  flooding  of  the  GDL  is  one  of  the  most  investigated  phe¬ 
nomena  in  numerical  work  but  Yamada  et  al.  [91]  performed 
experimental  analysis  for  the  extent  of  water  flooding  of  the  cath¬ 
ode  GDL  with  switching  interdigitated  and  conventional  flow  field 
designs  to  measure  pressure  drop  and  concluded  that  water  flood¬ 
ing  is  a  direct  function  of  the  wetting  properties  of  the  catalyst  layer 
and  the  GDL.  Additionally  Benziger  et  al.  [92]  measured  the  resis¬ 
tance  to  the  flow  of  water  through  the  GDL  by  applying  hydrostatic 
pressure  across  the  GDL  and  the  effects  of  pressure  on  the  water 
transport  through  the  voids  were  analyzed  in  detail.  It  was  con¬ 
cluded  that  water  flowed  through  only  1%  of  the  void  volume  in 
the  GDL.  Lin  and  Nguyen  [93]  also  employed  pressure  drop  calcu¬ 
lations  to  measure  the  effect  of  GDL  thickness  and  hydrophobicity 
on  flooding  levels.  An  optimum  condition  was  suggested  for  both 
water  flooding  levels  and  the  species  transport  through  the  porous 
GDL  because  the  hydrophobic  pores  support  the  gaseous  transport 
while  the  hydrophilicity  aided  the  liquid  water  transport.  It  was  also 
noticed  that  employing  micro-porous  layer  (MPL)  in-between  the 
catalyst  layer  and  GDL  reduced  the  water  flooding  which  makes  it 
possible  to  use  thinner  GDLs  for  same  operation  performance  level 
of  PEM  fuel  cell.  Litster  et  al.  [94]  developed  a  novel  approach  to 
visually  study  the  water  flooding  levels  in  the  porous  GDL  by  using 
fluorescence  microscopy  technique.  The  movement  of  water  was 
monitored  by  following  the  light  emitting  dye  with  a  microscope 
fitted  with  CCD  camera.  It  was  shown  that  liquid  water  is  trans¬ 
ported  by  fingering  and  channeling  similar  to  the  blotting  paper 
effect. 

In  many  CFD  analysis  of  PEM  fuel  cells,  the  flooding  of  the  gas 
channels  is  usually  ignored  keeping  in  view  the  fact  that  liquid 
portion  makes  up  only  minute  fraction  of  the  total  volume  of  the 
gas  channel.  However,  it  is  very  desirable  for  liquid  water  in  gas 
channels  to  be  removed  because  of  two  adverse  effects  that  can 
severely  degrade  the  effective  operation  of  PEM  fuel  cells,  i.e.,  flood¬ 
ing  in  channels  cause  a  liquid  film  to  cover  the  GDL  surface  and 
hinders  in  optimum  transport  of  gaseous  species  inside  GDL.  And, 
removal  of  liquid  water  at  shutdown  may  not  be  performed  ade¬ 
quately  [88].  Kumbur  et  al.  [95]  have  experimentally  investigated 
the  droplet  behavior  and  instability  in  rectangular  PEM  fuel  cell 
channels  because  the  physics  of  detachment  size  and  behavior  of 
water  droplets  plays  an  effective  role  in  the  liquid  water  removal 
from  a  surface.  The  findings  of  this  work  indicated  a  strong  relation 
between  contact  angles  (wettability)  and  the  departure  diameter 
on  each  other  along  with  the  surface  roughness  of  the  GDL  and  the 
characteristics  of  the  fluid  flow  through  the  channel.  Zhang  et  al. 
[60]  have  also  performed  similar  investigations  but  mainly  focused 
on  the  water  accumulations  at  the  geometric  corners  of  the  flow 
field  at  low  flow  rates. 

4.3.  Case  study  of  liquid  water  simulation  in  a  cathode  of  PEM 
fuel  cell 

A  two  dimensional  model  of  a  cathode  of  PEM  fuel  cell  has  been 
simulated  to  study  the  effect  of  liquid  water  on  the  performance 
in  terms  of  current  density  and  voltage.  The  reason  of  selecting 
the  cathode  side  only  is  based  on  the  fact  that  the  oxygen  reduc¬ 
tion  reaction  (ORR)  is  the  rate  limiting  reaction  in  PEM  fuel  cell 
electro-chemistry  [65].  The  salient  features  of  this  work  are  that 


7910 


M.A.  Khan  et  al.  /  Journal  of  Power  Sources  196(2011)  7899-7916 


Fig.  3.  Contours  of  water  saturation  at  different  current  loads;  (a)  0.2  A  cm-2,  (b)  0.5  A  cm-2,  (c)  0.8  A  cm-2  (d)  the  calculation  domain. 


agglomerate  catalyst  model  was  used  to  predict  the  reaction  rate 
and  the  effect  of  Knudsen  diffusion  was  also  included  alongside 
multi-component  diffusion  of  gaseous  species  as  given  in  Eqs.  (8) 
and  (9).  To  include  the  effect  of  temperature,  LTNE  approach  was 
utilized  because  it  has  been  previously  established  [46,48]  that 
fluid  temperature  is  lower  than  predicted  by  one  equation  model 
because  of  higher  thermal  conductivity  of  the  solid  structure  of  the 
media.  All  the  simulations  were  performed  on  Fluent®  with  3rd 
order  of  spatial  dicretization  of  the  domain.  In  this  case,  both  the 
inlet  and  the  out  of  the  domain  were  simulated  as  infinite  sink  for 
liquid  water,  i.e.,  zero  saturation  value  at  the  both  the  boundaries 
of  the  working  domain  (Fig.  3d). 

The  simulated  results  of  the  water  saturation  for  three  load  lev¬ 
els  i.e.,  0.2,  0.5  and  0.8  (A cm-2)  are  presented  in  Fig.  3(a,  b,  and 
c).  In  all  cases  it  can  be  seen  that  despite  of  zero  water  saturation 
values  at  the  inlet  and  outlet,  the  catalyst  layer  is  the  most  affected 
part  having  higher  levels  of  water  saturation  as  compared  to  GDL 
that  can  cause  significant  loss  in  efficiency  of  overall  process  by 
making  reaction  sites  unavailable  for  the  reactions.  Also,  it  can  be 
seen  that  by  increasing  the  current  density  or  load  the  saturation 
level  also  increases  suggesting  higher  reaction  rates  and  water  pro¬ 
duction  according  to  Eq.  (2).  Other  sources  that  contribute  to  the 
overall  generation  of  water  are  condensation  of  water  vapors  and 
the  electro-osmotic  drag  that  can  move  one  to  five  water  molecules 
from  anode  to  cathode  per  ion  transfer  [96].  Although,  operation  at 
low  loads  inherently  safe  guards  the  cathode  from  water  flooding 
but  at  such  conditions  a  significant  loss  in  protonic  conductivity 
was  also  observed  as  compared  to  high  load  operation  of  PEM  fuel 
cell. 

To  validate  the  model,  the  results  of  this  case  study  were  com¬ 
pared  to  those  produced  by  Sun  et  al.  [97]  in  which  the  agglomerate 
model  for  electro-chemical  reactions  was  developed  to  study  the 
influence  of  the  structural  parameters  on  the  catalyst  layer.  As  seen 


in  Fig.  4  that  both  the  predictions  almost  coincide  at  low  current 
densities  but  as  the  current  density  is  increased  a  gradual  deviation 
is  observed  between  the  two  cases  where  the  model  presented  by 
Sun  et  al.  over-predicts  the  current  density.  This  deviation  between 
the  two  cases  can  be  attributed  to  the  effects  of  water  saturation 
because  the  model  presented  by  Sun  et  al.  was  single  phase  only  i.e., 
no  liquid  water  was  assumed  to  exist  in  the  computational  domain. 
So,  it  can  be  concluded  that  one  of  the  effects  of  liquid  water  inside 
calculation  domain  is  to  limit  the  reaction  kinetics  by  effectively 
covering  the  reaction  sites  and  blocking  the  path  of  the  reactant 
gaseous  flow. 


Fig.  4.  Comparison  of  single  and  two-phase  models  for  PEM  fuel  cells. 


M.A  Khan  et  al.  /  Journal  of  Power  Sources  196(2011)  7899-7916 


7911 


4.4.  Two-phase  model  discrepancies 

The  list  of  two-phase  models  for  PEM  fuel  cells  can  continue  to 
an  infinite  number.  Many  authors,  as  reviewed  before,  have  pro¬ 
duced  novel  models  and  validated  their  results  against  previously 
established  or  experimental  work.  With  increase  in  the  compu¬ 
tational  power,  the  numerical  modeling  of  the  PEM  fuel  cells  has 
outpaced  experimental  work.  However,  experimental  work  is  still 
considered  the  ultimate  test  for  all  CFD  models  and  for  thorough 
validation  it  is  necessary  that  a  comparison  check  should  be  made 
between  theoretical  (CFD)  and  experimental  findings.  Keeping  this 
in  view,  Mench  [16]  has  recently  compared  the  CFD  work  with 
high-resolution  neutron  and  X-ray  imaging  results  of  liquid  water 
distribution  by  Turhan  et  al.  [98],  Weber  and  Hickner  [99]  and 
Hartnig  et  al.  [100]  and  found  out  that  various  assumptions  and 
theories  about  liquid  water  distribution  are  actually  quite  miss- 
leading  under  normal  operating  conditions  and  development  of 
new  models  was  highly  urged.  The  miss-match  of  experimental 
and  theoretical  findings  can  be  attributed  to  the  reasons  discussed 
below  (for  detailed  discussion  the  readers  is  referred  to  [16],  this 
section  only  summarizes  the  latest  findings  in  the  same  reference). 

4.4.1.  Large  liquid  water  accumulation  at  CL/MPL  interface 

Two  phenomena  have  been  pointed  out  by  Mench  [16]  that 
result  in  accumulation  of  liquid  water  in  large  quantities  at  the 
CL/MPL  interface; 

i)  Due  to  the  surface  roughness  of  mating  catalyst  layer  and  MPL, 
there  is  always  a  chance  that  large  voids  regions  are  left  behind 
even  at  high  commercial  scale  compactness. 

ii)  Cracks  found  at  the  stated  interface  are  highly  probable.  As 
reported  by  Mench  [16],  cracks  of  20  [xm  wide  may  run  through 
MPL  and  CL  and  constitute  total  of  8%  of  the  surface  areas. 

In  CFD  modeling  the  above  two  anomalies  are  generally 
neglected  and  a  perfect  CL/MPL  interface  is  assumed.  Recent  X-ray 
images  [100]  have  shown  a  large  quantity  of  water  accumulation 
at  the  interface  and  sometimes  they  represent  5-20%  of  the  total 
liquid  water.  Furthermore,  the  presence  of  large  scale  cracks  alters 
the  flow  pattern  of  liquid  water  by  providing  large  flows  through 
these  cracks. 


4.4.3.  Bruggeman’s  relation 

Usually  the  effect  of  saturation  in  the  porous  media  is  approx¬ 
imated  by  Bruggeman’s  relation  that  only  restricts  the  flow  of 
species  in  terms  of  reduced  porosity,  as; 

Deff  =  D,x  [eT(l— s)]  (11) 

According  to  Eq.  (1 1 ),  the  only  effect  by  the  presence  of  liquid 
water  is  to  reduce  diffusivity  of  the  gaseous  species  by  reducing  the 
void  volume  in  porous  media.  However,  as  noted  by  Mench  [16], 
this  relation  ignores  the  effect  of  blocking  of  pathways  for  gaseous 
flow  through  porous  media  that  can  significantly  affect  the  fuel  cell 
performance. 


4.4.4.  Other  discrepancies  in  two-phase  models 

Apart  from  the  previously  explained  ignorance  reigning  the 
two-phase  models,  few  more  considerations  need  to  be  taken  into 
account  for  vivid  description  and  picture  of  saturation  in  PEM  fuel 
cells,  as  outlined  here; 


i)  The  inter-phase  change  between  water  humidified  gaseous 
species  and  liquid  water  is  usually  modeled  as  [8,23]; 
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The  evaporation/condensation  process  in  PEM  fuel  cells  is  gov¬ 
erned  by  the  saturation  pressure  of  the  water  vapor  and  the 
temperature.  The  validity  of  this  process  is  still  a  question  mark 
as  it  has  not  been  experimentally  verified  and  experimental 
works  suggest  that  even  fractional  difference  in  temperature  can 
have  significant  deviation  from  the  simulated  results  by  Eq.  (12) 
[16]. 

ii)  Modeling  of  liquid  water  distribution,  usually  a  porous  jump 
condition  is  assumed  at  the  interfaces  of  components  because 
of  the  difference  in  porosity  and  radius  of  pores  in  two  adja¬ 
cent  materials  to  account  for  the  continuity.  These  assumptions 
predict  that  GDL  has  more  capacity  to  store  liquid  water  as  com¬ 
pared  to  MPL  because  of  its  higher  porosity,  but,  experimental 
results  show  the  other  way,  i.e.,  GDL  represents  low  accumula¬ 
tion  region  for  liquid  water  [16]. 


4.4.2.  Applicability  of  Lev  ere  tt  function  for  liquid  water 
distribution 

Not  all,  but  many  two-phase  models  proposed  for  water  man¬ 
agement  in  PEM  fuel  cells  utilize  the  Leverett  function  that  gives 
dominance  to  capillary  flow  in  the  porous  media,  given  as  [101  ]; 


where  Pc,  cr ,  s  and  I<  represent  the  capillary  pressure,  interfacial 
tension  between  phases,  porosity  and  the  absolute  permeability  of 
the  medium,  respectively.  However,  the  proposed  Leverett  func¬ 
tion  (also  called  J-function)  was  initially  proposed  for  flow  of  water 
through  soil,  whereas,  the  porous  media  used  for  PEM  fuel  cell 
construction  is  highly  compact  and  anisotropic  in  nature.  Kumbur 
et  al.  [102]  has  shown  that  significant  deviation  in  direct  measure¬ 
ment  of  water  quantity  as  compared  to  the  proposed  distribution 
by  the  Leverett  function.  So,  it  is  very  urgent  that  new  and  well 
descriptive  governing  equations  are  developed  for  water  flow  in 
PEM  fuel  cell  porous  media.  Recently,  Kumbur  et  al.  has  proposed 
new  relationships  for  governing  the  capillary  flow  in  PEM  porous 
media  [103-105],  based  on  direct  measurement  by  experimental 
techniques. 


As  discussed  above,  despite  of  active  model  developments  for 
PEM  fuel  cells  and  understanding  the  water  management  issue, 
many  inter-related  phenomena  still  need  to  be  accounted  and  a 
vigorous  and  robust  approach  is  required  that  can  encircle  all  the 
PEM  fuel  cell  issues  and  produce  efficient  and  stable  results  that 
can  help  in  commercialization  and  help  PEM  fuel  cells  to  compete 
with  highly  developed  conventional  counterparts  that  have  taken 
centuries  to  grow  and  groom. 

5.  Meso-scale  models  for  the  cathode 

The  electrodes  of  a  PEM  fuel  cell  are  volumetric  and  are  com¬ 
posed  of  catalyst  layer  and  GDL.  One  of  the  key  components  of 
PEM  fuel  cell  is  the  catalyst  layer  where  electro-chemical  reac¬ 
tions  take  place.  The  catalyst  layer  is  often  made  of  porous  mixture 
of  carbon  supported  platinum  particles  (Pt|C),  ion  and  electron 
conducting  material  (Nation),  and  also  provides  the  transport  of 
reactant  gaseous  species  to  the  reaction  sites.  Each  electrode  in 
PEM  fuel  cells  is  separated  by  an  ion  conducting  electrolyte  that  is 
usually  20-100  |xm  thick  [106].  The  Pt|C  particles  have  the  typical 
dimensions  of  lOOnm  and  are  thoroughly  distributed  in  the  cat¬ 
alyst  layer.  Furthermore,  Pt|C  particles  are  covered  with  the  same 
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ionic  conducting  material  (~10-7m  thick  [107])  as  used  in  elec¬ 
trolyte  preparation  to  insure  three-phase  contact.  Additionally,  for 
even  distribution  of  the  gaseous  species  and  the  electron  conduc¬ 
tion  from/to  the  catalyst  layer,  a  GDL  is  incorporated  to  a  thickness 
of  mm.  On  component  and  system  scale,  the  dimensions  are  ranged 
between  10  and  100  cm  and  more  [108].  In  order  to  encompass 
all  the  length  scales  in  PEM  fuel  cells,  a  multi-scale  model  needs 
to  be  developed  that  has  much  shorter  range  (beyond  continuum 
approximation)  than  current  CFD  scales  that  are  mostly  based  on 
Navier-Stokes.  For  categorization  of  the  fluid  flow  regime,  the  most 
commonly  tool  used  is  the  Knudsen  number  defined  as; 

fen  =  j  (13) 

where  X  is  the  mean  free  path  of  the  molecules,  and  l  represents 
the  characteristic  flow  dimension.  The  selection  of  the  character¬ 
istic  dimension  depends  on  the  length  scales  of  the  gradients  of 
pressure,  density,  velocity  and  temperature.  The  Knudsen  number 
is  useful  for  determining  whether  statistical  mechanics  or  the  con¬ 
tinuum  mechanics  formulation  of  fluid  dynamics  should  be  valid: 
if  the  Knudsen  number  is  near  or  greater  than  one,  the  mean  free 
path  of  a  molecule  is  comparable  to  or  larger  than  the  length  scale 
of  the  problem,  and  the  continuum  assumption  of  fluid  mechanics 
is  no  longer  a  good  approximation.  In  recent  years,  for  PEM  fuel  cell 
multi-scale  modeling  many  authors  have  opted  the  Lattice  Boltz¬ 
mann  (LB)  method  which  is  applicable  over  the  whole  range  of 
Knudsen  number,  including  the  continuum  regime.  The  LB  method 
provides  a  better  alternative  to  conventional  CFD  analysis  of  fluid 
flow  for  deeper  insight. 

5  A.  Lattice  Boltzmann  method  for  fluid  flows 

In  recent  years,  Lattice  Boltzmann  (LB)  method  has  emerged  as 
an  alternative  simulation  tool  for  predicting  fluid  flows  and  has 
been  found  very  successful  for  interfacial  dynamics  and  complex 
boundaries.  Lattice  Boltzmann  method  is  based  on  microscopic 
models  and  mesoscopic  kinetic  equations,  unlike  conventional 
numerical  methods  in  which  macroscopic  equations  are  discretized 
over  the  spatial  domain  [109].  The  underlying  theory  of  the  Lattice 
Boltzmann  method  is  to  develop  simplified  kinetic  models  cap¬ 
turing  key  features  of  micro-  or  mesoscopic  physics  so  that  the 
desired  macroscopic  equations  are  satisfied  since  macroscopic  fluid 
flow  is  a  collection  of  many  microscopic  particles  in  the  system 
[109,110].  Here,  the  LB  method  is  briefly  presented  for  basic  know¬ 
how,  for  detailed  understanding  one  can  refer  to  Succi  [111],  Chen 
and  Doolen  [109],  Sukop  and  Throne  [112],  Park  et  al.  [113],  and 
Spaid  and  Phelan  [114]. 

The  Lattice  Boltzmann  method  is  a  mesoscopic  method  that  is 
considered  in-between  the  continuum  based  technique  and  the 
molecular  dynamics  technique  which  handles  individual  particles 
in  the  flow  field.  In  this  method,  a  large  group  of  molecules  are 
assumed  to  move  about  a  lattice  and  collide  with  each  other.  At 
each  lattice  point,  the  particles  translate  at  discrete  directions  and 
all  particles  in  one  direction  are  grouped  together.  Generally  LB 
method  can  be  divided  into  two  sequential  steps  as  (i)  streaming 
and  (ii)  collision.  The  streaming  process  describes  the  movement 
of  a  particles  to  adjacent  lattice  point  in  the  direction  of  travel 
while  keeping  mass  and  momentum  as  conserved  quantities.  The 
intra  collision  of  particles  is  defined  by  the  collision  step  [115].  The 
general  governing  equation  for  the  LB  method  is  given  as  [109]; 

rij  (x+e,Z\x,  t  +  At)  =  n2-(x,  t)  +  (n(x,  t))  i  =  0,  1,  . . . ,  M  (14) 

The  Eq.  (14)  is  itself  quite  cumbersome  to  solve  because  of  an 
infinite  sets  of  particles  moving  along  different  directions  and  col¬ 
lisions  occurring  with  respect  to  the  scattering  rule.  A  modification 
is  generally  enforced  to  limit  only  one  particle  moving  in  certain 


direction  with  a  certain  velocity.  This  simplification  is  termed  as 
the  exclusion  principle  and  reduces  tracking  of  particles  to  a  finite 
and  manageable  number  for  a  given  time  step  [109].  For  scatter¬ 
ing  step,  which  is  non-linear  in  nature,  Higuera  and  Jimenez  [116] 
proposed  to  assume  the  distribution  close  to  equilibrium  state  to 
avoid  non-linearlization.  In  Eq.  (14),  nz  is  a  real  variable  represent¬ 
ing  the  mass  per  unit  volume  of  the  particle  translating  with  a  speed 
i,  while,  e2  is  the  direction  vector  and  M  represents  the  total  number 
of  lattice  points  in  consideration.  For  2D  geometries  M  ranges  from 
4  to  9,  whereas,  for  3D  cases  the  total  number  of  lattice  points  can 
be  either  15, 17  or  19. 

Since  1992,  a  considerable  attention  has  been  paid  to  use  LB 
methods  for  different  technologies.  According  to  the  data  collected 
by  Sukop  and  Throne  [112],  a  variety  of  fields  have  adopted  this 
technique  with  physics  and  mathematics  being  the  main  bearers. 
However,  only  a  few  papers  with  LB  method  and  PEM  fuel  cells 
as  title  or  keywords  have  been  produced  to  date.  Here  a  review  of 
these  is  presented  that  are  specifically  focused  on  PEM  fuel  cells. 

Generally  GDL  is  treated  as  homogeneous  and  isotropic  porous 
material  without  considering  the  manufacturing  details  of  the 
preparing  material.  Bundles  of  carbon  fibers  are  used  to  prepare  the 
carbon  papers  that  are  stacked  together  to  form  the  GDL.  Due  to  the 
structural  alignment  of  carbon  fibers  there  is  strong  dependency 
of  physical  properties  to  spatial  directions,  i.e.,  in-plane  (parallel  to 
surface)  and  through-plane  (perpendicular  to  surface).  As  described 
in  [  1 1 5  ],  GDL  can  be  aptly  assumed  to  consist  of  randomly  arranged 
cylinders  with  diameter  (7-12  pan)  being  very  small  as  compared 
to  their  length.  So,  due  to  inherent  alignment  formed  during  mak¬ 
ing  of  carbon  papers  for  PEM  fuel  cells  there  is  a  preferential  fiber 
direction.  Van  Doormaal  and  Pharoah  [115]  have  applied  LB  method 
to  determine  the  porosity,  permeability,  flow  direction  and  fiber 
directions  to  have  a  deeper  insight  of  GDL.  The  geometry  for  this 
study  was  generated  using  the  modified  Monte  Carlo  technique  by 
Himilton  [117].  The  alignment  of  fibers  was  varied  from  0  to  90°, 
i.e.,  from  parallel  to  perpendicular  directions  and  numbers  of  cases 
were  simulated  for  all  directions  with  porosity  ranging  from  0.6  to 
0.81 .  The  proposed  relation  between  porosity  and  the  permeability 
of  the  GDL  by  the  authors  is  given  as; 

s3.6 

0.26- - r2  in-plane 
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where  r  is  the  fiber  radius  and  predictions  were  found  to  be  within 
95%  of  the  prediction  by  Carmen-Kozeny  equation  [118]. 

Park  and  Li  [113,119]  have  also  performed  LB  simulations  for 
the  in-homogeneous  GDL  of  a  PEM  fuel  cell  with  considering 
multi-phase  flow  via  two  different  LB  formulations.  The  simulation 
domain  for  this  study  was  400  p,m2  consisting  of  voids  and  solid 
obstacles  represented  by  the  fiber  material.  It  was  demonstrated 
that  flow  squeezed  at  narrow  paths  in-between  the  fibers  causing 
pressure  drop  and  increase  in  velocity.  Additionally,  it  was  noted 
that  large  flow  blockages  can  be  caused  by  even  small  obstacles 
due  to  their  inherent  orientation  and  fibers  in  parallel  to  the  main 
flow  directions  had  no  significant  effect  on  flow  resistance  imply¬ 
ing  that  the  permeability  is  not  only  a  function  of  the  porosity  of 
the  GDL  material  but  it  is  greatly  affected  by  the  orientation  of  the 
fibers  also.  To  include  the  effect  of  liquid  water,  the  inter-particle 
interaction  model  was  also  used  as  LB  methods  are  capable  of  pro¬ 
viding  robust  predictions  of  the  interfaces  between  two  or  more 
phases  [119].  One  of  the  findings  by  the  authors  was  that  some  of 
the  liquid  water  moving  along  the  flow  direction  was  captured  in 
the  porous  media  and  it  can  be  concluded  that  the  captured  water 
can  cause  either  blocking  of  the  porous  path  or  covering  of  reac¬ 
tion  sites  if  it  is  in  the  catalyst  layer.  Similarly  it  was  also  predicted 
that  LB  method  is  far  superior  to  conventional  CFD  techniques 
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to  portray  unsteady  liquid  water  accumulation/removal  process. 
Karimi  and  Li  [120]  carried  out  numerical  calculations  for  a  single 
pore  of  the  membrane  material  to  investigate  the  electro-osmotic 
flow  through  Poisson-Boltzmann  and  Navier-Stokes  equations, 
and  predicted  that  the  pore  dimensions  have  significant  impact  on 
the  drag  coefficient  for  water  transport.  Similar  work  was  carried 
out  on  microstructure  of  the  membrane  material  by  Okada  et  al. 
[121]  to  demonstrate  the  dependence  of  the  ion  mobility  on  the 
microstructure  of  the  membrane  material. 

6.  Microscopic  models  for  the  catalyst  layers 

It  has  been  already  stated  in  the  previous  sections  that  the 
physio-chemical  features  of  the  catalyst  layer  is  one  of  the  per¬ 
formance  limiter  in  PEM  fuel  cells.  The  catalytic  activity  in  PEM 
fuel  cells  is  dependent  on  the  catalyst  particle  size,  shape,  distri¬ 
bution,  gaseous  access  to  the  reaction  sites,  ionic  and  electronic 
conductivity  and  thermal  distribution.  Catalyst  material,  process¬ 
ing  and  catalytic  reactions  are  all  inter-related  disciplines  and  a 
good  optimum  design  can  only  ensure  high  performance  catalyst 
layer  [122]. 

Liu  [122]  has  developed  top-down  and  bottom-up  approaches 
in  analysis  of  the  catalyst  materials.  In  top-down  approach,  a  multi¬ 
scale  analysis  approach  is  used  to  segregate  the  catalyst  material 
into  individual  technology  components  signifying  the  interaction 
of  each  to  the  overall  performance,  while  the  bottom-up  approach 
helps  in  developing  and  improving  the  catalyst  design  and  stability. 
Although,  with  advances  in  modern  science  and  technology,  higher 
degree  of  understanding  has  been  attained  for  the  catalyst  materi¬ 
als  and  performance,  but  still,  it  is  treated  as  a  black  box  in  most  of 
the  modeling  techniques.  A  need  is  felt  here  to  study  the  catalyst 
layer  at  particle  or  pore-size  levels. 


6.1.  Catalyst  layer  modeling  at  particle-  and  pore-size  level 

The  catalyst  layer  in  PEM  fuel  cells  forms  the  backbone  of  the 
unit  cells.  Apart  from  the  design  considerations  for  other  compo¬ 
nents,  the  optimization  of  the  catalyst  material  and  its  distribution 
inside  the  layer  is  very  essential  for  high  performance  and  utiliza¬ 
tion.  Recently  it  has  been  developed  that  manipulating  the  catalyst 
structures  can  significantly  enhance  its  performance  and  for  solid 
catalyst  particles  as  in  PEM  fuel  cells,  the  decrease  in  particle  size 
helps  in  significant  performance  upgrade  and  increasing  the  cata¬ 
lyst  size  can  cause  a  rapid  decline  in  performance  [122].  Since  Pt  is 
an  expensive  metal,  the  control  of  catalyst  design  is  very  eminent 
for  commercial  competition  of  PEM  fuel  cells  with  other  available 
energy  production  devices. 

Lee  and  Cho  [123]  determined,  through  chemical  quantum  cal¬ 
culations,  that  the  arrangement  of  individual  Pt  particle  can  lead 
to  the  optimization  of  catalyst  utility  and  performance.  Configura¬ 
tions  of  611  Pt  atoms  with  size  of  3.1  nm  have  been  found  to  be 
the  most  suitable.  However,  in  PEM  fuel  cells,  not  only  the  size  and 
configuration  are  important  but  transport  of  ions,  electrons  and 
reactants  to  the  catalytic  sites  is  also  very  essential.  In  PEM  fuel 
cells,  coating  of  larger  carbon  particles  with  nano-sized  Pt  parti¬ 
cles  submerged  in  electrolyte  material  has  been  termed  as  one  of 
the  most  effective  catalytic  distribution  patterns  for  high  utiliza¬ 
tion  [39,124].  The  parameter  used  to  describe  the  performance  of 
the  catalyst  in  PEM  fuel  cells  is  the  effective  surface  area  in  PEM  fuel 
cells  [125].  Using  finer  particle  size  leads  to  an  increased  effective 
surface  area  per  unit  volume  of  the  catalyst  particle  and  increasing 
of  the  catalyst  loading  (mass  per  unit  surface  area)  also  results  in 
increased  reaction  rates  [125]. 

Recently,  Siddique  and  Liu  [126]  have  digitally  constructed 
a  3D  model  of  the  catalyst  layer  at  nano-scale  using  controlled 


quasi-random  algorithms  numerically  duplicating  the  experimen¬ 
tal  fabrication  process  aiming  to  quantify  porous  nanostructure 
and  investigate  the  mechanism  of  nano-scale  electro-chemical 
reactions  and  percolation  networks.  It  was  concluded  that  a  thresh¬ 
old  percolation  exists  for  species  transport  through  the  catalyst 
structure  and  also  altering  the  optimum  value  of  the  agglomerate 
number  reduces  the  electro-chemical  reactions.  Similarly,  Lange 
et  al.  [127]  also  performed  3D  nano-scale  simulations  by  recon¬ 
structing  the  catalyst  layer  using  a  stochastic  approach.  The  focus 
in  this  work  was  to  compute  effective  transport  parameters  over 
selected  sets  of  operating  conditions  for  variety  of  microstructures. 
Additionally,  the  effects  of  water  vapor  and  temperature  profiles 
have  been  studied  in  depth.  For  the  validation  of  the  computed 
results,  a  detailed  comparison  has  also  been  performed  by  experi¬ 
mental  results  by  the  same  group  [128]. 


6.2.  Agglomerate  modeling  approach 


In  open  literature  many  mathematical  models  for  the  catalyst 
layer  have  been  proposed  from  zero  to  three  dimensional  models. 
Among  all,  the  flooded  agglomerate  model  is  the  most  descriptive 
and  predicting  model  [129].  In  this  model,  the  carbon  supported 
Pt  particles  in  form  of  agglomerate  are  immersed  in  a  thin  film  of 
electrolyte.  The  catalyst  layer  consists  of  a  network  of  intercon¬ 
nected  micro-  and  macro-sized  pores  through  which  the  gaseous 
species  reach  the  surface  of  the  agglomerates.  There  upon,  the  reac¬ 
tant  species  diffuse  through  a  thin  layer  of  the  electrolyte  to  reach 
the  reaction  site  [130].  The  agglomerate  model  has  been  able  to 
give  deeper  insight  into  the  physio-electro-chemical  phenomena 
in  simulations  and  modeling  of  the  catalyst  layer.  However,  the 
consideration  of  the  catalyst  layer  to  be  composed  of  carbon  sup¬ 
ported  Pt  with  flooded  electrolyte  film  as  a  continuum  medium  has 
made  it  difficult  to  analyze  the  discrete  distribution  of  the  catalyst 
phase  in  the  agglomerate  [131]. 

Yan  and  Wu  [  1 30],  Antonio  et  al.  [  1 32  ]  and  Bultel  et  al.  [133]  have 
developed  various  micro-models  for  the  catalyst  layers  in  which 
mass  and  ion  transfer  have  been  addressed  separately  by  treating 
the  agglomerate  and  the  electrolyte  material  (Nafion)  as  discrete 
and  segregated  components.  The  proposed  micro-models  have 
been  able  to  provide  deeper  insight  into  the  detailed  mass  and  ion 
transfer  mechanisms  at  pore  levels,  and  particle  size  relation  and 
dependence  in  the  overall  performance  for  the  electro-chemical 
reactions. 

The  relation  between  the  generation  of  the  current  per  unit 
volume  as  a  function  of  electric  and  ionic  potentials,  reactant  con¬ 
centrations  and  the  material  properties  of  the  catalyst  layer  is 
represented  by  a  modified  Butler-Volmer  kinetics  as  [134-136]; 


V-ii 


=  4  FC“k 


^film 


AaggDo2?nafion 


(16) 


where  S,  Aagg,  Dq2  nafion,  k  and  E  represent  the  thickness  of  Nafion 
covering  the  agglomerate,  agglomerate  area,  diffusivity  of  oxygen 
through  electrolyte,  reaction  rate  and  the  effectiveness  factor.  Con¬ 
sidering  the  first  order  reaction  kinetics,  the  analytical  expression 
for  the  effectiveness  factor  yielded  by  applying  the  mass  conserva¬ 
tion  equation  for  spherical  agglomerate  is  given  as  [136]; 

E  =  2^2  [30  cotli (30)  -  1  ]  (17) 

where  0  is  the  Theiles  modulus  for  the  system  and  is  expresses  by 

[8]; 


0  =  ? 


D, 


O2, nafion 


(18) 
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in  which,  f  is  the  characteristic  length  of  the  agglomerate  in  terms 
of  volume  per  unit  surface  area,  usually  given  as  Ragg/3  for  spher¬ 
ical  structures  as  used  in  preparation  of  catalyst  for  PEM  fuel 
cells  [136]. 

7.  Model  verification 

Since  all  numerical  studies  have  been  done  theoretically  by 
solving  a  combination  of  equations  backed  by  empirical  or  exper¬ 
imental  data  with  a  set  of  assumptions  to  make  it  feasible  with 
respect  to  both  time  and  computational  power  at  hand,  it  is  very 
important  to  validate  the  results  obtained.  Most  of  the  fuel  cells 
models  have  been  validated  against  the  measured  data  using  the 
polarization  curve  or  V-I  curve.  But  mere  comparison  with  V-I 
curve  does  not  ensure  complete  robustness  and  accuracy  of  the  pre¬ 
dicted  results  [  1 6].  According  to  Mench  [  1 6],  it  is  not  surprising  that 
a  good  agreement  between  the  model  and  a  simple  polarization 
curve  can  be  achieved  but  does  not  necessarily  validate  the  internal 
distribution  of  parameters  like  heat,  water,  and  charge,  since  under¬ 
prediction  of  one  parameter  may  be  balanced  by  over-prediction 
of  other  and  vice  versa.  Additionally,  fuel  cell  system  perfor¬ 
mance  fluctuates  with  different  electrode  assemblies,  however,  the 
numerical  simulations  are  not  flexible  enough  to  accommodate  all 
details  to  represent  the  fluctuations  in  the  performance  [16].  Sim¬ 
ilarly,  Pharoah  et  al.  [43]  has  also  reported  anomalies  in  model 
comparison  methods  using  the  V-I  curve  for  different  cases  while 
reviewing  the  material  anisotropic  effects  on  the  fuel  cell  perfor¬ 
mance.  It  was  noted  that  the  V-I  curve  for  both  anisotropic  and 
isotropic  electronic  conductivity  were  almost  identical,  however, 
the  internal  current  distributions  varied  significantly  for  different 
load  conditions.  For  the  same  load,  the  maximum  current  density 
occurred  under  the  rib  area  of  the  fuel  cell  for  the  isotropic  con¬ 
ditions  but  since  the  anisotropy  of  electrical  conductivity  makes 
the  conduction  much  higher  in  in-plane  direction,  it  was  noticed 
that  the  maximum  shifted  to  the  channel  center  line  in  the  latter 
case. 

Since  fuel  cell  predictions  are  dependent  on  a  number  of  param¬ 
eters,  and  most  of  which  have  not  yet  been  properly  characterized. 
Usually  to  match  the  predictions,  one  or  more  parameters  are  tuned 
to  obtain  matched  results  via  single  V-I  curve.  But,  proper  valida¬ 
tion  of  numerical  models  can  only  be  achieved  by  a  comparison 
with  detailed  and  local  experimental  data  and  results,  but  lack  of 
such  data  prohibits  the  proper  validation.  Since  fuel  cells  are  rela¬ 
tively  infant  in  age  comparatively,  there  is  a  tremendous  need  to 
develop  methods  and  techniques  for  real-time  data  acquisition  for 
validation  of  numerical  models.  However,  it  has  been  suggested 
by  Pharoah  et  al.  [43],  in  absence  of  real-time  data,  that  different 
V-I  curves  under  various  operating  conditions  should  be  com¬ 
pared.  Although  this  will  not  ensure  complete  validation  of  the 
models  but  operation  under  different  set  of  conditions  will  pos- 
turize  the  general  trend  in  the  performance  response  of  a  fuel  cell 
model. 

8.  Conclusions 

PEM  fuel  cells  represent  a  promising  future  for  the  energy  pro¬ 
duction  with  low  to  zero  greenhouse  gas  emissions.  Although, 
many  advancements  have  been  made  in  the  recent  years  in  PEM 
technology  both  in  terms  of  insight  into  internal  phenomena 
and  development,  however,  some  major  issues  still  need  to  be 
addressed  before  rendering  PEM  fuel  cells  for  large  scale  com¬ 
mercialization.  Among  many,  water  management  is  an  old  time 
problem  that  has  not  been  fully  understood  and  characterized. 
Since  PEM  fuel  cells  operate  under  different  load  conditions,  it  is 
quite  difficult  to  set  a  fixed  parameter  for  the  quantity  of  water 


as  at  higher  load  conditions  removal  is  deemed  necessary  to  avoid 
blocking  of  both  reaction  sites  and  pathways  for  gaseous  flow,  while 
in  low  operation  levels  less  water  quantity  exhibits  a  decreased 
ionic  conductivity  of  the  electrolyte.  So,  it  is  very  eminent  that  the 
water  management  issue  has  to  be  resolved  in  the  near  future  to 
make  PEM  fuel  cell  competent  with  respect  to  other  fuel  cell  types 
and  conventional  energy  sources.  With  the  advances  in  computa¬ 
tional  technologies  in  terms  of  both  speed  and  efficiency,  CFD  has 
become  an  optimum  tool  to  perform  detailed  study  of  PEM  fuel 
cells  in  all  aspects  from  understanding  the  internal  phenomena  to 
design  optimization.  A  great  deal  of  work  has  been  done  since  then 
with  a  variety  of  models  and  results.  Initially,  most  models  were 
limited  to  2D  and  single  phase  flows.  But,  with  deeper  interest 
in  managing  the  water  issue  inside  PEM,  later  models  developed 
were  extended  to  two-phase  with  both  2D  and  3D  geometries. 
But  still  a  complete  model  has  not  been  produced  yet  basically 
due  to  both  computational  limits  and  inherent  complex  phenom¬ 
ena  occurring  inside  PEM  fuel  cells.  The  discrepancies  in  modeling 
water  management  issue  still  needs  a  deeper  insight  ranging  from 
development  of  new  formulations  to  include  physical  and  chemical 
effects.  A  test  case  was  also  presented  to  study  the  effects  of  water 
saturation.  It  was  observed  that  the  catalyst  layer  was  the  most 
affected  area  as  maximum  water  saturation  effects  were  noted 
there  with  increasing  trends  with  increase  in  operating  current 
densities. 

Although  CFD  analysis  has  provided  an  overall  behavior  of  PEM 
fuel  cells  but  complete  validation  of  the  results  in  terms  of  robust¬ 
ness  and  capability  to  predict  the  transport  phenomena  still  needs 
to  be  verified  against  experimentally  established  results,  as  it  has 
been  often  observed  that  although  the  predicted  V-I  curves  have 
same  profiles  but  internal  distribution  of  parameters  like  current 
density,  heat  and  water  saturation  may  vary  significantly.  Also, 
PEM  fuel  cells  involve  multi-scale  parameters  and  phenomena,  e.g., 
the  electro-chemical  reactions  and  the  charge  transport  which  are 
best  evaluated  at  micro-  or  meso-scale  levels,  so,  the  multi-scale 
approach  will  further  elaborate  the  secretive  picture  of  the  fuel  cell 
operations.  Recent  developments  in  application  of  Lattice  Boltz¬ 
mann  method  to  flow  analysis  provides  the  facility  to  simulate  the 
flows  at  meso-scale,  although,  it  is  computationally  expensive  to 
apply  to  a  complete  single  cell  presently,  but  use  of  this  method  to 
PEM  fuel  cell  technology  is  pacing  up  and  some  researchers  have 
already  proven  its  worth  in  evaluating  physical  properties  such  as 
permeability  with  direction  dependent  profiles  when  applied  to 
selective  and  manageable  dimensions.  Further  development  in  this 
field  can  be  explored  to  reveal  more  insight  to  PEM  fuel  cells  and 
establish  robust  and  reliable  results. 
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